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ABSTRACT 

LeVeque and Yee recently investigated a one-dimensional scalar conservation law 
with stiff source terms modeling the reacting flow problems and discovered that for the 
very stiff case most of the current finite difference methods developed for non-reacting 
flows would produce wrong solutions when there is a propagating discontinuity. A 
numerical scheme, ENO/SRCI^ is proposed in this report for solving conservation laws 
with stiff source terms. This scheme is a modification of Harten’s ENO scheme with 
subcell resolution, ENO/SR. The locations of the discontinuities and the characteristic 
directions are essential in the design. Strang’s time-splitting method is used and time 
evolutions are done by advancing along the characteristics. Numerical experiment 
using this scheme shows excellent results on the model problem of LeVeque and Yee. 
Comparisons of the results of ENO, ENO/SR, and ENO/SRCD are also presented. 
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1. INTRODUCTION 


In the investigation of numerical methods for reacting flow problems, LeVeque 
and Yee [4] recently considered certain fundamental questions concerning the quality 
of numerical solutions. Namely, in extending current finite difference techniques de- 
veloped for non-reacting flows to reacting flows, can one: (i) develop stable methods, 
(ii) obtain “high resolution” results with sharp discontinuities and second order accu- 
racy in smooth regions, and (iii) obtain the correct jumps at the correct locations? 
They introduced and studied the following one-dimensional scalar conservation law 
with parameter-dependent source term 

Ut + U x = ( 1 ) 

^(u) = -n u (u - (u - 1), (2) 

where n is a parameter. This equation becomes stiff when the parameter fJ. is large. 
Although this linear advection equation with a source term represents only a simple 
model of reacting flow problems, by studying the numerical solutions one encounters 
some of the intriguing difficulties sure to occur in solving more realistic models. 

In their study, two different approaches were used to construct second order ac- 
curate numerical methods. One approach was to use a modification of MacCormack’s 
predictor-corrector method for conservation laws, together with two TVD-like versions 
with appropriate limiters. The other approach was based on the second order accurate 
Strang splitting method [5]. Their numerical tests revealed that stable and second 
order schemes can be devised by using either of these approaches. However, in study- 
ing the ability of these methods in dealing with propagating discontinuities, it was 
reported that for a fixed mesh and for the very stiff case, all the methods produced so- 
lutions that look reasonable and yet are completely wrong, because the discontinuities 
are in the wrong locations. Their investigation pointed out that the main difficulty is 
the smearing of the discontinuity in the spatial direction, which in turn introduced a 
nonequilibrium state into the calculation. To avoid this difficulty, it will be necessary 
to increase the resolution near the discontinuity, at least for the purpose of evaluating 
0(u). 

Integrating Eq. (l) over [xy_i,xy + i]x[i n ,t n +i], one obtains 



1 ft n+ 1 [X , 

+-r- / / if>(u(x,t))dxdt, (3) 

J tn J x ,_^ 

where u" denotes the cell average of u over [xy_j.,Xy + i] at t n . From this integral 
equation formulation, one can see that a source oJ error is also from the evaluation 
of the double integral term in (3). Because in most numerical methods the function 
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ip(u) is replaced by 0 evaluated at a fixed value of u. This may produce a reasonable 
approximation only when u is smooth. 

The purpose of this paper is to show that numerical methods can be devised 
to overcome the above mentioned difficulties. We will construct a numerical scheme 
which, when applied to Eq. (1), results in stable solutions with excellent resolutions 
at the correct locations of the discontinuities. Essential to the construction of this 
scheme is the application of Harten’s ENO reconstruction with subcell resolution [2]. 
The subcell resolution is an idea based on the observation that unlike point values, 
cell averages of a discontinuous piecewise smooth function contain information about 
the exact location of the discontinuity within the cell. Using this observation in his 
study of conservation laws, Harten [2] proposed reconstruction techniques and obtained 
modifications of the ENO schemes [3] showing significant improvement in the resolution 
of contact discontinuities. Basically, when good approximations to the exact locations 
of the discontinuities inside the cells can be obtained, it is then possible to have good 
reconstruction of the solution at each time step. Here we will also demonstrate that 
when the information on the location of the discontinuity is used in treating the source 
term, the results will improve significantly. 

To maintain second order accuracy, we will use the approach of Strang’s time- 
splitting method [5] in which one alternates between solving the conservation law with- 
out the source term and the ordinary differential equation modeling the chemistry. The 
numerical solution v n+1 at time step t n+1 is computed from v n by 

<-" +1 = S#(f)5/(Al)S*(^)v". (4) 

where Sf(At) represents the numerical solution operator for the conservation law with- 
out the source term over a time step At, and S^(^) represents the numerical solution 
operator for the ordinary differential equation 

u t = V'(u), (5) 

over At/2. In terms of the integral equation formulation (3), Sf takes care of the 
flux terms represented by the two single integrals at and S ^ handles the double 

integral term. 

We will outline the construction of the two operators Sf and 5^, in section 2, which 
depends on the characteristic directions as well as the ENO reconstruction with subcell 
resolution procedure [2]. For the purpose of comparison, we will also test numerical 
schemes which use a ENO and a ENO with subcell resolution as Sf. These algorithms 
will also be stated briefly in section 2. In section 3, we report the numerical results 
obtained from using the above schemes on the same model problem of LeVeque and 
Yee [4]. A conclusion will be given in section 4. 


2. CONSTRUCTION OF THE SCHEMES 

We first describe the structure of the operator Sf(At). At the time step t n , suppose 
that we have obtained the numerical solution v n = {t/”}, where represents an 
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approximation to the cell average u*. Then, to obtain S/(At) v n , we use the following 
steps: 

1. Obtain a reconstruction R(x\v n ) from the given values v n . 

2. Modify this reconstruction R(x;v n ) to obtain R{x\v n ) using the subcell resolution 
when discontinuity is detected. 

3. Advance R{x\ v n ) along the characteristics from t n to t n +i and then take cell 
averages to complete S/(Af) v n . 

In the scheme we propose, the steps 1 and 2 will follow the basic ENO reconstruc- 
tion procedure with subcell resolution of Harten [2]. The reconstructed solution func- 
tion R(x\ v n ) here is a piecewise quadratic polynomial obtained by using the primitive 
function approach. For the sake of completeness, we will describe in straightforward 
terras the procedures used. For more details and general discussions on reconstruction 
and subcell resolution, see Harten [2], 

Step 1. ENO Reconstruction 

Over each cell [xy_i,Zy + i], choose i = i(j) such that 

\ v i+2 - 2 <+l +v i I = rmn {| v jfc+2 - 2v k + 1 + V k \ : k = J- 2 ij ~ 1 >j}- 
Let Rj{x\v n ) denote the reconstructed quadratic polynomial over this cell. Then 


where 


Rj (*; v n ) = aj + s, {x - Xj) + - Cj (x - xjY 

c j = (v? +2 — + v ?)/(Ax) , 

s i = ( u ?+i - V i)/ Ax + U - * - \) c i Ax » 

ay = V* — Cj (Az) 2 /24. 


Step 2. Subcell Resolution 

To detect a discontinuity in a cell [xy_i,Xy + i], we define 


( 6 ) 


( 7 ) 


= ±\f' 


2 4 
i* I - 


Fj(z) = ~ [ / Rj-i(x-, v n )dx + J ’ * Rj+ i(x; v n ) dx) - Vy . 
*i- h * 

In the schemes we tested, the following criterion is used. If 

|sy| > |sy_i|, |sy | > |«y+i|, and Fy(x y _i) Fj{a r J + i) < 0, 
we consider that there is a discontinuity at Oj in this cell satisfying 

Fy(<?y)=0. 


( 8 ) 


( 9 ) 


(10) 
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The location Oj can be approximated by using any standard root-finding method. We 
simply use the bisection method in our experiment. 

Now, if there is a discontinuity inside the cell, a modified reconstruction I2y(x; v n ) 
is used , where 



Rj-i{x;v n ), 

Rj+i{x',v n ), 


x- i < X < Oj, 

J 2 

Oj <x<x j+h . 


Otherwise, we use 

Rj{x\v n ) = Rj(x-,v n ). 


Step 3. Time Evolution and Cell Averaging 

Here we choose to describe the scheme for the conservation law 


ut + au x — 0, a > 0. 


( 11 ) 


It can be handled similarly for a < 0. Consider the case a At < Ax and that there 
exists a discontinuity at Oj inside the cell [xy_i,xy + i] with 

Oj < Xj + 1 — a At, 

as shown in Fig. 1. The idea is that, following the characteristics from t n to t n+ \, we 
will find an approximation to 


S/(A!) v? = 



( 12 ) 


In our present scheme we use the following simple computation. Let x m and x p denote 
the midpoints in the intervals (xy i — a At,0j) and (0y,Xy + i — a At) respectively (see 
Fig. 1). Then compute 


S f {At) v if = [iEy_i(x m ; v n ) ( Oj - x y _i +aAf) 


+iEy + i(x p ;u n ) (x y+ i — a At — Oj)] /Ax. (13) 

Other locations of Oj and the cases with smooth regions can be treated similarly and 
easily. It is quite simple to modify the above scheme for more general equations and 
also to write higher order versions of it. 

Now, let us describe the operator Sy, (At). It is essentially the approximation of 
the double integral term in Eq.(3). Let us use the same notations introduced above 
and refer to Fig. 1. To advance the value I2y_x(x m ; t/ n ) from t n to t n+i , we again 
follow the characteristics to obtain an approximate value 


Rj—i(x m ',v ) + Atip(Rj-i(x m -,v )), 
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using the simple Euler’s method. Let z m and z p denote the midpoints in the intervals 
,$j) and (0y,Xy+i) respectively. Then, for the case a > 0 and Oj < x J+ i - a At, 
we use 

(*, - *,_$) 

+ ^(.Ry_i(z m ;t> n ) + AtV’(^y-i( a; m;w w ))) (*y - *y-$ + aAt) 

+ 0(fiy+i(* p ;v n )) (x y+ i - 0y) 

+ V»(.Ry+i(zp;t; n ) + A*V’Cffy+i(av»v n ))) fcy+i “ e i - a At)]. 


Again, other situations are handled similarly. 

The resulting algorithm then takes the following form: 

vr l = s *(f )SAAt)s ' f,i f )v ?- 


(14) 


(15) 


Two things are essential in the design of the above algorithm, namely, the location 
of the discontinuity Oj and the characteristic direction. We approximate Oj by following 
the ENO reconstruction with subcell resolution procedure. For this reason, we denote 
this algorithm by ENO/SRCD. In applying this algorithm to Eq. (1), one takes the 
value a = 1. In treating more general equations, the formulas in step 3 can be easily 
modified. 

For the purpose of comparison, we have tested several other schemes. We shall 
report the results from using a ENO scheme and also a ENO scheme with subcell 
resolution, ENO/SR, of Harten [2] as the operator 5/ (At). The following version of 
the ENO scheme has been used in [1]. 

ENO Scheme: 

For the operator Sj (At), we use 


where 


with 


/y+i — fj+i — 2 

+/ R (M z y+i>*n+iMy+i( I y+i>*»+i))]> 


(17) 


Az 


Vj{ x j+L,t n ) — v" + 2 5 y 


Ax 


®y+i(*y+j»*») - v y+i 2 5;+1 


Az 


vy(xy + i,t n +i) = Vj + ~~z~ *j - Atasj, 


Ax 


*V+i(*y+$ — w y+i — o s,+1 A f a 5 y+i> 
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where the sy’s used in the computation come from (7) in step 1, a = 1 for Eq.(l), and 
f R [vL,VR) denotes the flux at the origin in a Riemann problem with vl to the left and 
vr to the right. 

In [2], there is a simpler version of second order ENO scheme. The difference is 
from the approximation of the two single integrals at Zy ± i in Eq.(3). The ENO in [2] 
uses a midpoint rule to do the approximation, while a trapezoidal rule is used in the 
above formulation. Since in our test Eq.(17) produces slightly better computational 
results, we choose to report it here. The operator (At) will be the same as in Eq.(14) 
and the final algorithm also takes the same form as in Eq.(15). 

ENO/SR Scheme: 

The operator S/(At) is now replaced by Harten’s second order ENO scheme with 
subcell resolution [2]. The entire algorithm is denoted here also by ENO/SR. The 
construction of Sf(At) is described again for Eq.(ll) in the form of Eq.(16) with 


fj+ 


1 

3 


—ENO 

~ fj+k ^>+2 ’ 


where /y+i will be the same as in Eq.(17) and the correction term £y + 1 is computed 
as follows. If the discontinuity contion (9) is not satisfied, then 

9j + 1 = ^ (" ” !) ( 2u ~ !) c i ( Aa: ) 2 > v = a At/ Ax; 

( [(Ax -a At) (v" - o Atsy/2) - 6y_i(xy_|,xy + i - a At)] /At, 

_ I when Fj (x y + j - a At) Fj (x y _ ^ ) > 0, 

[6y + i(xy + i - a Af,Xy + i) - a At (v" + (Ax - o At) sy/2)]/Af, 

, otherwise, 

expression 6y(y 1 , 1 / 2 ) is used to mean 

rv a 

b]{yi,y 2 )= I Rj(x;v n )dx. 

Jyi 

In the above formulas, all the ay’s, sy’s, and cy’s come from (7) in step 1. We also use 
the same operator S 1 p(At) from (14) and the final algorithm again takes the formal 
form of (15). 


else 




and the 


3. COMPUTATIONAL RESULTS 


We use the same mesh and initial data as in the model problem of LeVeque and Yee 
[4] to test the ability of the above schemes in dealing with propagating discontinuities. 
Thus Eq.(l) is solved together with the initial condition 


u 



1, 

0, 


X < 0.3, 
x > 0.3. 
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if 

if 



We take Ax = 0.02, At = 0.015, and the domain in x to be from 0 to 1. For 
comprison with [4], we also show the results at t = 0.3 and for the cases ft = 1, 10, 
100, and 1000. Figure 2 shows the computed results using the. ENO, ENO/SR, and 
ENO/SRCD schemes for fi = 1, 10, and 100. For the very stiff case, f.i = 1000, both 
ENO and ENO/SR schemes fail to produce stable solutions for the above mesh in 
our computational experiment. Only the ENO/SRCD scheme still produces excellent 
results as shown in Fig. 3. However, when we reduce the size of At to one half of the 
original, i.e., At = 0.0075, and march 40 time steps, excellent results are again obtained 
from both ENO and ENO/SR schemes as shown in Fig. 4. Of course, reducing At 
means the reduction of the stiffness of the system. The difficulty arises from the fact 
that in both the ENO and ENO/SR schemes, the computation of the numerical flux 

—ENO * * . 

fj + L still produces “large ” error in the spatial direction. 

The computational results obtained here compare favorably to those in LeVeque 
and Yee [4]. 


4. CONCLUSIONS 

We have proposed a numerical scheme ENO/SRCD, which is a modification of 
Harten’s ENO scheme with subcell resolution ENO/SR, for solving conservation laws 
with stiff source terms. We use Strang’s time-splitting method and treat the conserva- 
tion law without the source term and an ordinary differential equation with the source 
term representing the chemistry sequentially. For both the conservation law solution 
operator Sj and the ordinary differential equation solution operator S^, the locations 
of the discontinuities and the characteristic directions are essential in their design. The 
main difference between the construction of this Sf and that of Harten’s ENO/SR is 
that the time evolution here is accomplished by advancing along the characteristics 

explicitly, while ENO/SR uses the numerical flux /y+i followed by a correction term 
based on the location of the discontinuity. The operator S ^ for the ordinary differen- 
tial equation also advances along the characteristics. Our numerical experiment using 
this scheme shows excellent results on the model problem of LeVeque and Yee [4] for 
reacting flows. Comparisons of the results of ENO, ENO/SR, and ENO/SRCD have 
also been presented. 
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ENO 


ENO/SR 


ENO/SRCD 



H = 10 



Fig. 2. Numerical results at t — 0.3 using ENO (first column), ENO/SR (second 
column), and ENO/SRCD (third column) schemes with discontinuous initial 
data for y = 1 (first row), y = 10 (second row), and y = 100 (third row). 
Ax = 0.02, At — 0.015, : true solution, ■ • •: computed solution 
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ENO/SRCD 
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Fig. 3. Numerical results at t = 0.3 using ENO/SRCD scheme with 
discontinuous initial data and \i = 1000. 

Ax = 0.02, At = 0.015, : true solution, • * computed solution 
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Fig. 4. Numerical results at t = 0.3 using ENO and ENO/SR schemes with 
discontinuous initial data and fi = 1000. 

Ax = 0.02, At = 0.0075, : true solution, • - •: computed solution 


12 



NASA 

National Aeronautics and 
Space Administralion 


Report Documentation Page 


1. Report No. 


NASA TM- 102384 
ICOMP-89-27 


2. Government Accession No. 


3. Recipient’s Catalog No. 


4. Title and Subtitle 

On the Application of Subcell Resolution to Conservation 
Laws With Stiff Source Terms 


5. Report Date 
November 1989 


6. Performing Organization Code 


7. Author(s) 

Shih-Hung Chang 


8. Performing Organization Report No. 

E-5123 


10. Work Unit No. 

505-62-21 


9. Performing Organization Name and Address 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135-3191 


11. Contract or Grant No 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546-0001 


13. Type of Report and Period Covered 
Technical Memorandum 


14. Sponsoring Agency Code 


15. Supplementary Notes 

Shih-Hung Chang, Department of Mathematics, Cleveland State University, Cleveland, Ohio 44115 and Institute 
for Computational Mechanics in Propulsion, Lewis Research Center, Cleveland, Ohio 44135 (work funded by 
Space Act Agreement C99066G). Space Act Monitor: Louis A, Povinelli. 


16. Abstract 

LeVeque and Yee recently investigated a one-dimensional scalar conservation law with stiff source terms modeling 
the reacting flow problems and discovered that for the very stiff case most of the current finite difference methods 
developed for non-reacting flows would produce wrong solutions when there is a propagating discontinuity. A 
numerical scheme, ENO/SRCD, is proposed in this report for solving conservation laws with stiff source terms. 
This scheme is a modification of Harten’s ENO scheme with subcell resolution, ENO/SR. The locations of the 
discontinuities and the characteristic directions are essential in the design. Strang’s time-splitting method is used 
and time evolutions are done by advancing along the characteristics. Numerical experiment using this scheme 
shows excellent results on the model problem of LeVeque and Yee. Comparisons of the results of ENO, ENO/SR, 
and ENO/SRCD are also presented. 


17. Key Words (Suggested by Author(s)) 

Subcell resolution; Conservation laws; End reconstruction; 
Characteristics direction; Stiff source team 


18. Distribution Statement 

Unclassified — Unlimited 
Subject Category 64 


19. Security Classif. (of this report) 

20. Security Classif. (of this page) 

21. No of pages 

22. Price* 

Unclassified 

Unclassified 

12 

A03 


NASA FORM 1626 OCT B6 


*For sale by the National Technical Information Service, Springfield, Virginia 22161 






